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Abstract 

We prove strong convergence of a semi-discrete finite difference method 
for the KdV and modified KdV equations. We extend existing results to 
non-smooth data (namely, in L^), without size restrictions. Our approach 
uses a fourth order (in space) stabilization term and a special conservative 
discretization of the nonlinear term. Convergence follows from a smooth- 
ing effect and energy estimates. We illustrate our results with numerical 
experiments, including a numerical investigation of an open problem re- 
lated to uniqueness posed by Y. Tsutsumi. 
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1 Introduction 

This paper is concerned with the study of a numerical approximation of the 
equation 

dtu + dlu + Pd^u^+^ ^Q, p^Q, fc = l,2, (1.1) 



with u = m(x, t), a; G M, t > 0. When fc = 1, the equation ( 1.1 ) is referred to as 
the Korteweg-de Vries (KdV) equation, and when fc = 2 as the modified KdV 
(mKdV) equation. 
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As is well-known, the KdV equation describes the unidirectional propagation 
of small-but-finite amplitude waves in a nonlinear dispersive medium. It ap- 
pears in several physical contexts, such as shallow water waves and ion-acoustic 
waves in a cold plasma. Also, the modified KdV equation has been used to 
describe acoustic waves and Alfen waves in plasmas without collisions. For a 
more complete description of the physical contexts concerning the Korteweg-de 
Vries equation and its generalizations, see [32] and the references therein. 

A large amount of work on the KdV equation was initially directed toward 
the study of solitary waves, i.e., solutions of the form u{x,t) = U{x — ct), 
especially the so-called soliton solutions, a class of solitary waves which preserve 
the form through nonlinear interaction (see [321 [T7] for surveys on solitons). One 
of the most relevant results in soliton theory was the development of the inverse 
scattering method, initially applied to the KdV equation by Gardner et. al. [7] 
and, in a general form, by Lax |14j . This technique was also used to obtain 
solutions of the KdV equation with low regularity [H [3 H] . 

Here, we concentrate on the numerical approximation of the solution of the 
Cauchy problem 



dtu + dlu + Pd^u''+^ =0, 13^0, fc = l,2, (1.2a) 
u{x,0) = ip{x), ipeL^. (1.2b) 



The mathematical problem of well-posedness for (1.2a),(1.2b) has been exten- 
sively studied. We refer to the pioneering results in [T1 [21121] and the improve- 
ments in [mill]. In these works, local well-posedness is proved in the Sobolev 
spaces 7?", s > 3/2, for generalized KdV (gKdV), in which the term d^u'^'^^ is 
replaced by dxV{u). 

Existence and uniqueness was also obtained in jS] [S] with initial data in 
weighed and spaces. In our numerical approach, we follow the energy 
method used in those papers. 

More recently, following the introduction by Bourgain [3] of certain Fourier 
spaces, the well-posedness result is strongly improved for data in negative Sobolev 
spaces (see the monograph [T^ and the references therein), and uniqueness of 
solution in is proved in 

Regarding the numerical solution of the KdV equation, convergence results 
have been proven for a linearized equation [TU] and for smooth solutions [TO] . 
However (to our knowledge) , the problem of proving rigorous convergence of nu- 
merical schemes without smoothness assumptions has only attracted attention 
in more recent years. Nixon [TH] proves the convergence of approximate solu- 
tions for a discretized version of gKdV, but for small initial data, only. That 
work is the numerical counterpart of [I3] . Finally we refer to the recent work by 
Pazoto et. al. ^D], dealing with the numerical treatment of the mKdV equation 
with critical exponent and a damping term, which shares some techniques with 
the present work. 

Thus, to the best of our knowledge, the problem of rigorous convergence of 
numerical schemes for the KdV and mKdV equations with general data in 
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has remained unsolved. The purpose of this paper is to fill that gap. 

Although the techniques we use to prove our convergence result are based 
on the ones in [9] (namely, the use of a fourth order stabilization term), their 
application to the numerical case is not trivial. Indeed, it is essential to use 
a special non-conservative discretization of the nonlinear term in (1.1). This 
idea dates back at least to HUj, and is also used in [2U]. Moreover, to obtain 
the necessary estimates for the numerical approximation, additional technical 
difficulties related to interpolators are encountered, with which we deal below. 

An outline of the paper follows. After some notations and definitions in 
Section [2] we prove our main convergence result in Section [3] In Section |4] we 
present some numerical experiments to illustrate our convergence results and 
test the accuracy of our scheme. 

Finally, in Section [Sj we investigate numerically an open question posed by 
Y. Tsutsumi [23] relating to the uniqueness of solution to the Cauchy problem 
for the KdV equation with measure initial data. This is done by means of the 
Miura transformation (see [13]), which relates solutions of the KdV equation 
with measure initial data to solutions of the mKdV equation with LP' initial 
data. As explained in more detail in Section [5] the numerical evidence we 
provide suggests that the Cauchy problem for the KdV equation with measure 
initial data is ill-posed. Note that, importantly, these numerical simulations 
involve discontinuous initial data in only, and, as such, are not covered by 
previous convergence results. 



2 Notations and definitions 

Let h denote our discretization parameter. We denote by = Uj{t) the (semi- 
discrete) difference approximation of u{xj,t), Xj — jh, j € Z. For /i > 0, we 
define the Banach spaces 

£1{Z) = {(z,) : z, e C, ^ ^\'^\' < 

For p = 2, we denote the usual scalar product by 

{z,w)h = y^^hzjWj, 

z = (zj), w — (wj). Let us also introduce the following standard notations for 
finite difference operators. For u — (uj), 

AhUj = D+D^Uj = D^D+Uj = — (uj+i — 2uj + uj^i), 
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Also, denote the translation operators by 

We obtain the following formulas for the discrete differentiation of a product, 

D+{vu) = vD+u + u+D+v (2.1a) 

D_{vu) ^ vD_u + u_D^v (2.1b) 

Do{vu) = vDqu + ^{u+D+v + u-D^v) (2.1c) 

Ah{vu) = vAhU + uAhV + D+vD+u + D^vD^u. (2.1d) 
Also, the difference operators verify in £f^ 

{D+u, v)h = -{u, D_v)h, {DqU, v)h = -{u, Dov)h, 

and so 

{AhU,v)h = (u, Ahv)h- 
We will also need to denote for a sequence (uj) and for a function w 

\]h\<R 

\\w\\p.R = \\w\\lp{-r,r), < R< +00. 
Finally, we introduce the continuous piecewise linear interpolator 

P^^u{x) = Uj + {x- Xj)'^^^ —, X e {xj,Xj+i), (2.2) 

Xj+i — Xj 

and the piecewise constant interpolator {Pqu){x) =Uj, a; G {xj,Xj-^i). 

3 Convergence results 

In this section, we prove our main result, Theorem |3.3[ which establishes the 



convergence of a numerical approximation of problem (1.2) 



Let us consider the semi-discrete finite difference scheme 

^u'^ + D^u'^ + p^^^liu'^fD^u'' + Do(v!'f+^\ + hAlu'' = 0, (3.1a) 
at k + 2 

u^(0) = (3.1b) 

where D^uj = D+D^D^Uj and A\ — D+D-D+D- denotes the difference bi- 
laplacian, and denotes the unknown grid function (u'j)j^z, '"j (i) being the 



approximation of the solution of (1.2 1 at the point {xj,t) 
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The term hAf^u^ is introduced in our scheme in order to obtain the uniform 
(in h) stabihty estimates necessary for the convergence proof. This term corre- 
sponds to the parabohc regularization used in [8] for the continuous problem. 
Also, the formally consistent discretization 



(3.2) 



is based on a corresponding one in \lQi and is also essential in our proof. See 
also [50] for an application of the same idea in a different setting. 
The following first result holds. 

Proposition 3.1. Let h > 0. Then, for each initial data ip^ £ ^h(^)' there 
exists a unique global solution u^{t) G C(IR;i?^(Z)) of (3.1 1. 



Proof. Existence of a unique local solution in follows from the Banach 

fixed-point theorem. The global existence is an immediate consequence of the 

□ 



uniform bounds on the R norm, established below in Lemma 



3.5 



Let Pi denote the continuous piecewise linear interpolator and let (p G 



be the initial data for the problem ( 1.2 1. Also, we denote by Cw{T, X) the space 



of weakly continuous functions from the interval I to the Banach space X. 
It is now convenient to introduce the notion of weak solution to the Cauchy 



problem ( 1.2a), ( 1.2b ). 



Definition 3.2. Let ip E L^(IR). A function u{x,t), a; £ M, i > 0, is a solution 
to the Cauchy problem ( 1.2a),( |L2b ) if 

1. ueL^,iiO,^);LHR)), 



2. For each test function (p E 2)(M x (0,oo)) one has 



u{dt4> + die/)) + pu''+'d^(f)dxdt = 0, 



(3.3) 



3. u{t) ipin L2(M) as t a.e. 

We now state the main result of this paper. 

Theorem 3.3 (Convergence of approximate solutions). Let p^ E he the 

initial data for the discretized problem (3.1), such that P^ip^ (p in L^(M) 
when /i — > 0. Then, for each T > 0, the sequence P^u^ satisfies 

P^u^^u m L'^{[Q,T];Hl^{^)) weak, 
P^u'^^u in L^{[0,T];Ll^ 



with u a solution of (1.2). Moreover, u satisfies 

«e (L°°nC„)([0,T];LLW), 
||"(i)||2 < ||V5||2 and ||u(i) — (p||2 — > as /i — ^ 0. 
The proof will be postponed to the next section. 



(3.4) 



5 



3.1 Main estimates 

First, let us record some inequalities which will be of use throughout. Let 

V = (vj) € £1{Z), q e [l,oo]. From we derive 

WPi'v - Pl^vWq^R < Ch\\Dnv\\q^R^h, < i? < OO , 

\\Pi'v~P^v\\g^R<C\\v\\qM,h, 0<i?<c5o, ge[l,+cx)], 
for some C independent of h. As a consequence, we obtain 

\\Piv\\g,R < C\\P^'v\\g..R = C\\v\\g,Rj,,yv e £UZ), < i? < +00. 
We will need the following inequalities. 



Lemma 3.4. Let 6 £ il(Z). Then, 



< CU\0\D^<I>\\1(', 



II0IIOO < 2ll^±<^lli.'^ 



(3.5) 
(3.6) 

(3.7) 
(3.8) 



Proof. The inequality (3.7 1 is a consequence of the Gagliardo-Niremberg in- 
equaUty and = PS^+- 

Moo = llPiVlloo < c||PiV||^/^||vPiV||^/^ 



while (3.8) is a consequence of the (continuous) inequality ||0||oo < ^ll^'lli- D 
We are now ready to state our first stability estimate. 



Lemma 3.5. Let u^{t) be a solution of (3.1 1. Then, for all t > 0, it holds 

\\u''{t)h.h < y%,h. (3.9) 



In particular, u is a global solution of (3.1 1. 

Proof. In the next proofs, we will, for simplicity, omit h from the notation. Take 
the £^-scalar product of the equation (3. la I with u = it'' to get 



Now, 



and so 



+ h{Alu,u) = 0. 

{Alu,u) = {AhU,Ahu) = \\Ahu\\l 
{v!'Dou,u) = {u^+^,Dou) = -{Dov!'+^,u), 

d 



dV 



L{t)\\i = -h\\Ahu\\i<0, 



(3.10) 



from which the conclusion follows. Notice how the discretization (3.2 1 leads to 
the non-increase of the norm. □ 
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The next lemma is a fundamental identity which, as we will see in Proposi- 
tion 3.7 below, implies a smoothing effect inherent to the equation: even though 
the initial data is only in ^^(Z), the solution of (3.1 ) is actually in a more regular 
space (uniformly in /i), namely, H^^^. 

Let p : M — )■ R be a bounded, strictly increasing, smooth function, with all 
its derivatives bounded. Write pj — p{xj), j G Z. For simplicity, we do not 
distinguish in our notation the continuous and the discrete p. 



Lemma 3.6. Let = u^{t) he the solution of the discrete problem (3.1 1. Then, 
u^ satisfies the identity 

^|b'/'«'ll2 + {D^u\D+pD^u^) + \{D+u\D^pD+u^) 



h{D+D^u'\pD+D^u'' 



'-{D+D^u'\D+pD_u'') 



+ -{D+u^D^p, D+D^u'^) - {D^u'',u'lDoD^p) 

~ h{D+D^u'\ D^pD^u'^) - h{D+D^u'\ D+pD+u'^) 
- h{D+D^u^,u''D+D^p) 

+ ^l^{in''r'^<D+p + ulD^p). 



(3.11) 



Proof. We take (3.1a), multiply by hpjUj, and sum over j e Z to obtain 



^ '^:\\p'/'u\\l + iD'u,pu)+(5^^[iu''Dou,pu)+{Dou''+\pu)]+hiAlu,pu) = 0. 



2dt' 



(3.12) 



We find from (2.1 ) 

{D^u,pu) = {D+DaD^u,pu) = -{DqD^u, D^{pu)) 



'{DoD_u,pD_u) - {DoD_u, u_D_p) ^ A + B. 



Since 



A = {D^u,pDqD^u) + {D^u, -{D+p{D^u)+ + D^p{D^u)^)), 



we obtain 



and so 



2A ^ {D^u, ^D+pD+u) + {D^u, ^D^p{D^u)^) 



- [{D^u, D+pD+u) + {D+u, D+pD^u)] 



{D+u, D+pD^u) 



= - [{D+u - D^u, D+pD^u) + {D^u, D+pD^u)] 
= ^{D+D^u,D+pD^u) + ]^{D^u,D+pD^u). 
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Similarly, 



= (D^u.u^DoD^p) + ^{D^u,D+pD^u) + ^{D+u,D^uD^p). 

Since 

{D+u,D^pD^u) = {D+u,D^pD+u) - h{D+u, D^pD+D^u), 
we obtain 

B = (D_u,u_DoD_p) + -{D_u,D+pD_u) 
+ ^{D+u,D_pD+u) - ^{D+u,D_pD+D_u). 



For the term in (3.12) corresponding to the discrete bi-laplacian, we derive 

{Alu,pu) = -{D^D+D^u,p^D^u + uD^p) 

= -{D^D+D^u,p^D^u) - {D^D+D^u,uD^p) 
= {D+D_u,pD+D_u) + {D+D_u,D_pD_u) 
+ {D+D^u, D+pD+u) + {D+D^u, uD+D^p). 

As to the remaining term in ( |3.12[ ), we find 

[u^'Dou^pu) + {Dqu''+\pu) - {Dau,pu''+^) - {v^+\Dq{pu)) 



k+l 



'■ 2 



{u+D+p + u^D^p)) . 



All these results together give (3.11). This completes the proof of Lemma 3.6 



□ 

As a consequence of the two preceding lemmas, we now prove the following 
result which states that, at the discrete level, u'^ is in H^^^. 



Proposition 3.7. Letu^ be solution of the discretized problem (3.1) with initial 
data (p = (p^ ^ t\{'L). Then, for each T > and for each R > 0, there exists a 
constant C = C{R,T, ||v||2,;i) such that, for all h> 0, 



I h\D±u''j\^ dt<c. 

Jo , , , , ^ „ 



(3.13) 



|j7i|<_R 



Proof. We apply Lemma |3.6| with a bounded, strictly increasing, smooth func- 
tion p, with all its derivatives bounded, and such that, moreover, p{x) > 1 for 
all X, and p'{x) = 1 for x G [—R,R]. Let us rewrite the identity (3.11), with 
obvious notation, as 

Ify^'^'Wl + {D-u\D+pD^u'^) + \{D+u\D^pD+u'^) 
+ h{D+D_u'',pD+D_u'') ^ Ai + --- + At. 
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Now observe that under our assumptions on p, the terms on the left-hand side 
(except the first) are non-negative, so we must bound the terms Ai. 
The terms Ai and A2 are similar and yield 

'^\{D+D^u'\D+pD^u'^)\ < '^\\{D+D^u'^)p^/%\\p-^/^D+pD^u% 

for all 77 > 0, where we have used the properties of p. Also, the terms A4 and 
A5 are similar and give 

hiD+D^u'\D^pD^u^) 

= h{D+D^u'\ D^pD+u'') - h^{D+D_u'\ D.pD+D.u'') = Bi + B2. 

The term Bi is similar to Ai, while 

IB2I < Ch^(D+D_u'\pD+D_u''). 

For the term ^3, we remark that 

{D^u'\u'lDoD^p) = {D^u'' - D^u'l,u'lDoD^p) + {D^u'!i,u'lDaD^p) 
= h{D^D^u'\u'LDoD^p) - {u'^,D^u''DoD^p) - {u'l,u''D+DoD^p) 

and so 

\{D_u\utDoD^p)\ < Ch\\p^^^D+D^u%\\u% + C\\uY2 
< Cih^\\p'/^D+D^u''\\l + C2\\u''\\l 

The term Aq = h{D^D^u''^,u'''D+D-p) is easily estimated using the Cauchy- 
Schwarz inequality. Finally, for the last term 

Ar^^l^^{{u'r^\ulD^p + u'^D^p), 

let us consider only k — 2, since the case /c = 1 is easier. Setting q = {D+pY^^ 
and u ^ , we obtain 

\{u\u+D+p)\<C\\u\\l\\qu^^ 
and from ( |3.8[ ) it follows 

llgu^lloo < \\{D^q)u^_ + q{uD^u + u^D^u)\\i < C\\u\\l + !|9i:>_it||2||M||2. 
Hence, 

\iu\u+D+p)\<C\\u\\i + C\\u\\l\\{D+py^'D^uh 

<e\\{D+py/^D^u\\l + -\\u\\t + C\\u\\i 
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for any e > 0. 

Choosing e, rj small enough, and for some small enough ho (which does not 
depend on R,T or (p), we find, after integrating (3.14) on [0,r] and using the 
previous estimates, 



±u'n'dt<C{R,T,M2). 



\]h\<R 



This completes the proof of Proposition 3.7 



□ 



3.2 Proof of Theorem Q 



The proof of Theorem |3.3| relies on Aubin's compactness result, which we state 
here, in a simplified form, for the reader's convenience. 

Lemma 3.8 f |151 p. 58]). Let \ < p < oo, T > 0, and consider reflexive Banach 
spaces Bq C B C Bi such that Bq is compactly embedded in B. Then, the space 

[v:v€L'{0,T-Bo), ^ G LP(0, T; } 

is compactly embedded in (O, T; _B) . 

In order to apply Lemma |3.8t we will use the following estimates. 



Lemma 3.9. Let be given by (3. la), (3. lb) and let T,R > 0. Then, there 
exists p > 1 such that 



h„,h\\'^ 



\P>' 



dt < C, 



(3.15) 
(3.16) 



uniformly in h, with C = C{R,T, ||(p||2.?i). 



Proof. First of all, note that since dxPi = PqD^, it follows from Lemma 3.5 



and Proposition 3.7 that the estimate (3.15) holds for each R> 0. 

Let us now prove the estimate (3.16). Let u'* be given by (3. la), (3. lb). We 
apply the piecewise linear continuous interpolator Pi to the equation (3.1a) to 
obtain 



P^'u" + P^D^u^' 
k + I 



dt 



(3.17) 



with — u^{t). We begin by estimating the term PiD^u'\ for which it is 
convenient to consider the decomposition Pi ~ {Pi — Pq) + Pq and analyze 
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the two resulting terms. For each test function G T){—R,R) we have 
\jh\<R '■" 



rxj+i 

y D^DoD^Uj / (j){x) dx 

h\<R. ^1 

22 D~Uj / T^W^ - 2^1) - 0(a; -h)- (j){x) + (t>{x + h)) dx 



\jh\<R 

and so (by Taylor expansion of (p) 

|(Po''i?V,0)| <C h\D^Uj\U"\\^ 

bh\<R 

<c[Y. /^l^-^jf) UU-{~B.,Ry 

Hence, by Proposition |3.7[ we have 

''||P^W||^_3(_^,^,)dt<C. (3.18) 
Next, if a; G {xj, a^j+i) and Vj = D^uj, we easily find 

{Pl^-P„')v'{x) = {x-x,f-^±^, 
and so, with obvious notation, 

((Pi'' - P^)D^v!', D+D^uj / {x- Xj)<p{x) dt 

= ^{D+D^u,,A,) 

A straightforward computation gives 



DoD_D_Aj 



{(j){x + h)~ 2(j){x) + 2(j){x - 2h) - (j){x - 3h)){x - x-j) dx 



2h^ 



from which we obtain by Taylor expansion of (j) and Proposition 3.7 

i-T 



I \{{P^-P',')D^u\<j>)\<c[ ^\^-''o\W"\Udt 



\jh\<R 



Ujl^l dt 



< C\\(t)\\HH-RM)- 
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From this and (3.18) we obtain the estimate 







|P^Z?V||^-.(_^_^) dt<C. 



In an entirely similar way, we arrive at 



Jo 



(3.19) 



(3.20) 



It remains to estimate the nonlinear terms in (3.17). Choose a smooth 
function 6* : M M such that e{x) = 1 if |a;| < i? and Oix) = if |a;| > i? + 1. 
Using (3.5) and (3.71 we derive, for i? > 0, fc = 1,2, 

\K[{^")'D„u'%^^^^ ,^ < \\Dou%^n+iM\\0A\kh 

<||i?o^i''||2.fl+i,/.(C + Cp+«"||^g;l/^) 

<C + C\\DW'&Xl^l 

with C — C(||(y3||2,/i, i?). Choosing p = 12/{3k + 5) > 1, we obtain from Propo- 
sition |321 



2,R+l,h — 



^ \\P^'[{u'')''Dou''']\\l^^dt<C + C \\D+U 



Since IKP/' - Fo'')"''ll3/2 < C\\Pl'u%/2 (cf. we conclude that 

£ \\Pi[iu'')''Dou^]\\l^^dt < C, fc = 1,2. 



(3.21) 



For the remaining nonlinear term PiDq{u^)*'^^^ we split it as above into 
(P/' - Po) + Pa - First, note that 

(3.22) 

and, by (3.6 1, the same estimate is obtained for ||Pi''(u'')'^^^||2,j?,- Next, since 
Pl^Do = d^Pl', we obtain for fc = 1, 2 



\piiDo{u^r^'\\'^'. 



dt 



< 



\d,pHu'^)'+^i;\_^j,^dt 

\PHu'r^'\\Tndt 



<C + C r \\D+uX,R+i,hdt<C. 
Jo 
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We now need to estimate (P^ — Pq)Do{u^)'^^^ . With computations similar to 
the ones after (3.181, we find 



\jh\<B, 

<C\\P,\u'^)'^ + \j<j>\\H^-R,R^, 



and by ( |3.22[ ), 







\{Pi,^ - pi^)Do{u'^rtH'^ 



dt<C + C 



\D+u'^\\lj,^^,^,dt<C. 



Thus we conclude that 







1,2. 



(3.23) 



The desired esti mate (|3.16 ), with p = 12/(3fc + 5) > 1, now follows from the 
estimat es (|3.17| , ( |3.19 ), ( |3.20] ), ( |3.21| , and ( |3.23p . This completes the proof of 

□ 



Lemma 13.91 

Proof of Theorem \3.3{ In view of the estimates in Lemma |3.9[ we apply Lem- 



3.8 withp = 12/(3fc + 5), Bo = H^{-R,R), B = Li{~R,R), q £ (l,oo], and 



Bi — H^^ {--R, R) (note that Bq C B with compact embedding). We conclude 
that, up to a subsequence, P^u''' converges weakly in L^{0,T; H^{~R, R)) and 
strongly in L?{Q,T] L'^{—R, R)). Using a diagonal argument, we obtain for a 
further subsequence 



P^u^^u in L2([0,r];ffi(-i?,i?)) weak. 



(3.24) 



h„.h 



P'^U 



in L2([0,r];L'?(-i?,i?)), i? > 0, ge(l,oo] 



for some u e L°°(0, T; ^^(k)) n L'^{Q,T]Hl^{W)), as /i ^ 0. Also, from ( [3^ 1 
and Proposition |3 . 7| we can conclude that 



(3.25) 



P^u'^^u in L^{0,T;Li{-R,R)), l<q<2. 



Now we must prove that u is a weak solution of the problem ( 1.2a),( |L2b ), 
in the sense of Definition |3.2| Let us apply the piecewise constant interpolator 
Pq to the discrete equation (3.1a|: 



hPl^Alu'' = 0. 



(3.26) 

First, consider the linear terms. We take a test function </> G ©(M x (0, oo)) and 
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compute in the sense of distributions 



(/){x, t) dx dt 



{dl<t) + 0{K)) dx dt 



= -{p^u'\dU) + 0{h) ^ -{u,dU) = {dlu,4>) 

as h ^ Q. The term hP^Af^u'^ is treated similarly and tends to zero as ft — 
in the sense of distributions. 

Now consider the nonlinear terms. Note that P^'Dq = d^P^ and write 
^0' = {Pa - Pi) + Pi - Using ^ and ^ we find 

||(P^ - < Ch\\D+iu^)'+%,i,,^ < Ch{l + \\D+u,\\l('^^,^,) 

and so 



ips-pniu') 



li\/,,h\k + l 







(3.27) 



in L'*/3(0,T;L^^M)) as /i -J> 0. Since Pq commutes with the nonlinearity, it 
follows from (3.25) that 

^ u'^+i in L\0,T;Ll^{R)). (3.28) 

Hence, we deduce from (3. 27), ( [3^ I that 

P/^Dq{u'')''+^ = d^P^iu'')''+^ ^ 

in the sense of distributions. For the other nonlinear term, we note that 
P^[{u'')''Dou^] = P^{u'')''P^Dqu''. We have 

P^iu'^f-^u'^ in L\0,T;Ll,{R)) 



and, from (3.24), 



Therefore, 



P^Dou'' = d^Pl'u'' - in i2(o, T; LL(K))- 



Po''[(^/'^)^-i^o^''] in ii(0,T;Li„,(M)). 



Multiplying (3.26) by a test function in D{R x (0, 00)), the above convergences 



allow us to conclude that u verifies the property (3.3) of Definition 3.2 

It remains to prove the weak L^(M)-valued continuity property, (3.4), and 
that u{t) — > in L^(M) as t — > a.e. To prove the weak continuity property, 
we remark that, for e I»(M), t e [0, T), 



{P,^u\t + p)~P^u\t),<j^)^ ( 



~df'' 



ds 
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and, from (3.16), we get 

\{Pi'u'^{t + p)-P^u\t),^)\<Cp 

with C ~ C (T, (p) , and so the family t (t) is miiformly bounded in L 

(see (3.6)) and weakly equicontinuous. Therefore, the Ascoli-Arzela Theorem 
implies that u € C„([0, cx)); L2(K)). 

Finally, since ||u(t)||2 < \\'p\\2 a-S- in t, we have 

\\uit)-^\\l<\\uit)\\l^2{uit),ip) + Ml 
<2Ml~2{uit),ip) ^0 



for almost every i — > 0+. This completes the proof of Theorem 3.3 



□ 



4 Numerical experiments 

4.1 A fully discrete, fully implicit scheme 

In this section, we present some numerical experiments to test the accuracy of 
our scheme and to illustrate our results. In order to fully discretize the semi- 



discrete equations (3.1a), we use a fully implicit Euler scheme, as follows. Given 



a time step r and a space step h, solve for each n — 1,2, . . . the equations 
,,"+1 _ i. 1 1 

+ D\r' + j^.ii-rrDour' + D^iu^^r^^] ^^^^^ 

+ 7]hAlu';+^ =0. 

We have set f3 — 1 and introduced a new viscosity parameter ?/ > allowing us 
to explicitly control the amount of viscosity in the scheme. 

As is standard in the numerical simulation of dispersive equations, we con- 
sider a sufficiently large spatial domain and initial data exponentially small 
outside some bounded region, ensuring that spurious wave reflection at the 
boundary of the domain remains negligible. 



Written in full, the scheme (4.1) reads 



,,"+1 _ 1 



+ 2P - + ^^r' - ^-l-l + "}^) = 0. 

Due to the nonlinear terms, it is necessary to perform a Newton iteration at 
each time step, which we carry out with a tolerance of 10^^ in the simulations 
below. To solve the pentadiagonal linear system at each iteration of Newton's 
method, we employ a standard LU decomposition method. 
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Relative L»2 error: || u_exact|L2 / 1| u - u_exact|L2, T = 10, tau = 0.001, eta = 1, k = 1. 















Relat 


ve error, 


k=i -g- 






































































3 E 


3 E 



500 1000 1500 2000 2500 3000 3500 4000 4500 5000 

Number of points 

Figure 1: Relative error as a function of the number of spatial points, k=l 
(KdV equation). 



4.2 Comparison with exact solutions 

The first step is to test our scheme with the known soliton solutions of the KdV 



equation (1.1) [H p. 140]. These read 



u{x,t) = [^{k + 2)chech'{'^c{x - ch))]^^'" 



(4.2) 



for arbitrary c > and consist of traveling waves with speed c^. We observe in 



passing that these exact solutions actually solve the equation (1.1 1 and not the 
slightly different version in [TBI P- 139]. 

In Figures [l] and [2j we present the error between the exact solution ( |4.2[ ) 
and the computed solution, at T = 10, computed on the domain x € (10, 50), as 
a function of the number of spatial points, for different values of the time step 
T, and, respectively, for fc = 1 and k — 2. 

One advantage of the present method is that it allows direct control of the 
amount of dissipation my means of the parameter rj in (4.1 ). We first note that 
our convergence results remain valid for any 77 > (but not for rj — 0). As 
would be expected, reducing the value of i] provides a sharper, less dissipative 
approximation. This is confirmed by our simulations, and in Figure[3]we present 
the error at T = 10 for various values of rj. Interestingly, setting 77 = 
sometimes provides a very good approximation, but not always, which is perhaps 
a consequence of the instability of the scheme without dissipation. 
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Relative L"2 error: || u_exact||_2 / 1| u - u_exact||_2, T = 10, tau = 0.001, eta = 1, k = 2. 
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Figure 2: Relative error as a function of the number of spatial points, k=2 
(mKdV equation). 
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Relative L»2 error: || u_exact||_2 / 1| u - u_exact||_2, T = 10, tau = 0.0001, i< = 1. Logscaie in y. 
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Figure 3: Relative error as a function of the viscosity parameter r] (KdV 
equation). 
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5 On an open question of Y. Tsutsumi 



In [23], the Cauchy problem for the KdV equation (1.1) with measure initial 
data is considered. In that work, the author addresses the open question of 
uniqueness of solution to the Cauchy problem for the KdV equation with mea- 
sure initial data in the following way. 

It is well known that a solution of the mKdV equation with initial data 
may be transformed, by the Miura transform u i— > M(u) — dxU + w^, into a 
solution of the KdV equation with a measure as initial data. Now, the family 
of functions 

c+1 

ul(.)^{ + ^ (5.1) 

^— , a;<0, 
X + c 

with c < —1, all verify M(u2) — (5(0), where 5 denotes the Dirac delta. There- 
fore, if Uc(x, t) is the solution of the mKdV equation with initial data vPc.{^), the 
question arises whether the Miura transform maps each of these different solu- 
tions to the same solution of the KdV equation with 5{Q) as initial data, or if, on 
the contrary, M(wc(a;,i)) varies with c, which would establish non-uniqueness. 
If the latter case is observed numerically, it would support the conjecture that 
the Cauchy problem for the KdV equation with measure initial data does not 
enjoy the uniqueness property. 

We have investigated this question numerically, with a high degree of preci- 
sion, and found that our numerical experiments support this lack of uniqueness 
conjecture. Thus, we have considered the mKdV equation with initial data given 
by M° (see (5.1)) for various values of c < —1, computed the solution Uc{x,t) 
up to some time T > 0, applied the Miura transform JA{uc{x,t)), and finally 
compared the solutions obtained. 

We have observed a clear dependence of M(uc(a;, i)) as c < —1 varies, see Fig- 
ure |4] This provides strong numerical evidence in support of a non-uniqueness 
property for the KdV equation with measure initial data and also a non-trivial 
test of the robustness of our numerical method: recall that the initial data (5.1 ) 
are discontinuous functions in only. 

Note that these simulations were computed with an accuracy of 30000 spatial 
points. Due to the slow decay of the solution, the computational domain is taken 
to be the interval [—500, 500], giving a value of /i = 1/30. We have also performed 
the computations with a coarser grid of 5000 points, and have found that the 
(natural) slight variation with h does not affect the overall qualitative aspect 
of the solution. In other words, the lack of uniqueness conjecture is strongly 
supported by our precise numerical experiments. 

Finally, we have verified as well that the result does not depend on the viscos- 
ity parameter rj appearing in (4.1). The simulations presented take rj = 0.001, 
but considering larger values of r/ (up to = 0.1) gives virtually indistinguish- 
able results. 
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Figure 4: Miura transform of solution for various values of c. T — 5, 30000 
spatial points, r = 0.001. 



In fact, it is easy to check that the more general family 

-, a; > 0, 



^Y^^^^^' (5.2) 

, X < 

a; + c 

verifies M(u° ^) = N{c, e)<5o, with N — {c + e — 1) / c. The same remarks about 
uniqueness apply, and so as a last test we have carried out simulations with 
(c, e) = (-1/4,1/4) and (c, e) = (-1,-2) (for which iV(c, e) = 4), performing 
the same comparison of the Miura transform of the computed solutions. 

For these simulations we have taken a very fine grid of 50000 spatial points, 
which corresponds toh = 0.02. The viscosity parameter is -q — 0.001. In Figure[5] 
we plot the Miura transform of the solution for two different values of (c, e), with 
T = 10, and 10000 and 50000 spatial points for each value of (c, e). 

Again, some variation with h is observed for the same values of (c, e), which 
is natural since the scheme includes dissipation. But the main thing to note 
are the appearance of two distinct solutions, one for each set of values of the 
pair (c, e), clearly apparent in Figure [5j The same distinction between the two 
solutions is also apparent for intermediate values of the number of grid points, 
whose solutions are seen to lie smoothly between the ones presented here. 

We can therefore conclude that our numerical experiments strongly indicate 
lack of uniqueness for the Cauchy problem for the KdV equation with a measure 
initial data. 
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Figure 5: Miura transform of solution for various values of c, e. T = 10, 10000 
and 50000 spatial points, r = 0.001. 
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